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Abstract 

We develop two classes of composite moment-free numerical quadratures for computing highly oscil¬ 
latory integrals having integrablc singularities and stationary points. The first class of the quadrature 
rules has a polynomial order of convergence and the second class has an exponential order of convergence. 

We first modify the moment-free Filon-type method for the oscillatory integrals without a singularity or 
a stationary point to accelerate their convergence. We then consider the oscillatory integrals without a 
singularity or a stationary point and then those with singularities and stationary points. The compos¬ 
ite quadrature rules are developed based on partitioning the integration domain according to the wave 
number and the singularity of the integrand. The integral defined on a subinterval has either a weak 
singularity without rapid oscillation or oscillation without a singularity. The classical quadrature rules 
for weakly singular integrals using graded points are employed for the singular integral without rapid 
oscillation and the modified moment-free Filon-type method is used for the oscillatory integrals without 
a singularity. Unlike the existing methods, the proposed methods do not have to compute the inverse of 
the oscillator. Numerical experiments are presented to demonstrate the approximation accuracy and the 
computational efficiency of the proposed methods. Numerical results show that the proposed methods 
outperform methods published most recently. 

Key Words: oscillatory integrals, algebraic singularities, stationary points, moment-free Filon-type 
method, graded points. 

1 Introduction 

We consider in this paper numerical computation of highly oscillatory integrals defined on a bounded in¬ 
terval whose integrands have the form /e 1K9 , where the wave number k is large, the amplitude function / 
may have weak singularities, and the oscillator g has stationary points of certain order. Computing highly 
oscillatory integrals is of importance in wide application areas ranging from quantum chemistry, computer¬ 
ized tomography, electrodynamics and fluid mechanics. For a large wave number k, the integrands oscillate 
rapidly and cancel themselves over most of the range. Cancelation dose not occur in the neighborhoods of 
critical points of the integrand (the endpoints of the integration domain and the stationary points of the 
oscillator). Efficiency of a quadrature of highly oscillatory integrals depends on the behavior of functions 
/ and g near the critical points. Traditional methods for evaluating oscillatory integrals become expensive 
when the wave number k is large, since the number of the evaluations of the integrand used grows linearly 
with the wave number k in order to obtain certain order of accuracy. The calculation of the integrals is 
widely perceived as a challenge issue. Calculating oscillatory integrals requires special effort. 

The interest in the highly oscillatory integrals has led to much progress in developing numerical quadra¬ 
ture formulas for computing these integrals. In the literature, there are mainly four classes of methods for 
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the computation: asymptotic methods, Filon-type methods, Levin-type methods and numerical steepest 
descent methods. The basis for convergence analysis of these quadrature rules is the asymptotic expansion 
of the oscillatory integral, an asymptotic expansion in negative powers of the wave number k. The leading 
terms in the asymptotic expansion may be derived from integration by parts m for the case when the 
oscillator has no stationary point. For the case when the oscillator has stationary points, the main tool is 
the method of stationary phase [241 [Ml . For a fixed wave number, the convergence order of the asymptotic 
method is rather low. To overcome this weakness, the Filon-type methods imunusni were proposed, which 
replace the amplitude function / by a suitable interpolating function. In many situations the convergence 
order of the Filon-type methods is significantly higher than that of the asymptotic methods. A thorough 
qualitative understanding of these methods and the analysis of their convergence order may be found in 
mmm for the univariate case and in m for the multivariate case. In these methods, interpolation 
at the Chebyshev points ensures convergence [23]. A drawback of the Filon-type methods is that they re¬ 
quire to compute the moments, which themselves are oscillatory integrals. For the cases having nonlinear 
oscillators, it is not always possible to compute the moments exactly. In [6] [7] f27i [37] , the moment-free 
Filon-type methods were developed. An entirely different approach without computing the moments is 
the Levin collocation method [21] . The Levin-type methods [22], 26] [28] [29 j reduce computation of the 
oscillatory integral to a simple problem of finding the antiderivative F of the integrand, where F satis¬ 
fies the differential equation F' + ing'F = f. The Filon-type methods and the Levin-type methods with 
polynomials bases are identical for the cases having the linear oscillator but not for the cases having the 
nonlinear oscillator [2U EH |38]. Numerical steepest descent methods [2] El [El for removing the oscillation 
converts the real integration interval into a path in the complex plane, with a standard quadrature method 
used to calculate the resulting complex integral. 

Although many methods were proposed for computing oscillatory integrals in the literatures, there are 
still a big room for improving their approximation accuracy and computational efficiency. The Filon-type 
and Levin-type methods require interpolating the derivatives of the amplitude / at critical points in order 
to achieve a higher convergence order. Even though computing derivatives can be avoided by allowing the 
interpolation points to approach the critical points as the wave number increases for the formula proposed 
in m, the moments cannot always be explicitly computed. In particular, certain special functions were 
used for calculating the oscillatory integrals in the case when / has singularities and g has stationary 
points. The formulas proposed in mm do not require computing the special functions and the moments of 
these formulas can be computed exactly, at the expenses of computing the inverse of the oscillator, which 
takes up much computing time. Numerical steepest descent methods also require computing the inverse of 
the oscillator or high order derivatives of the integrand [2j. 

The purpose of this paper is to develop efficient composite quadrature rules for computing highly 
oscillatory integrals with singularities and stationary points. The methods to be described in this paper 
require neither computing the moments of the integrand, inverting the oscillator, nor calculating the 
derivatives of / and those of g 1 . The main idea used here is to divide the integration interval into subintervals 
according to the wave number k and the singularity of the integrand. To avoid using the special functions, 
we first split the integration interval into the subintervals according to the singularity of / and the stationary 
points of g such that the integrand on the subintervals either has a weak singularity but no oscillation, 
or has oscillation but no singularity or stationary point. The weakly singular integrals are calculated by 
the classical quadratures using graded points [H H9U35]. To avoid using the derivatives of / and those of 
g' and avoid computing the inverse of the oscillator, we design a composite quadrature formulas using a 
partition of the subinterval formed according to the wave number k and the property of the oscillator g 
for the oscillatory integrals. These formulas can improve the approximation accuracy effectively, since the 
convergence order of the formulas computing the oscillatory integrals with smooth integrand and without 
stationary point may be increased by adding more internal interpolation nodes. Specifically, we develop 
two classes of composite moment-free quadrature formulas for the highly oscillatory integrals. Class one 
uses a fixed number of quadrature nodes in each subinterval and has a polynomial order of convergence. 
This class of formulas are stable and easy to implement. Class two uses variate numbers of quadrature 
nodes in the subintervals and achieves an exponential order of convergence. Convergence order of this class 
of formulas is higher than that of the first class. 


The quadrature formulas proposed in this paper have the following advantages. Comparing with the 
existing formulas, the proposed formulas need not computing the inverse of the nonlinear oscillator, or 
utilizing the incomplete Gamma function [1] for the oscillator with stationary points. These formulas not 
only reduce the computational complexity, but also enhance the approximation accuracy. The approxima¬ 
tion accuracy of these formulas is higher than that of the existing formulas for the case when the oscillator 
integral has stationary points and the oscillator is not easy to invert. 

We organize this paper in seven sections. In Section 2, we present an improved moment-free Filon- 
type method for the oscillatory integrals developed in m- In Section 3, we design a partition of the 
integration interval and propose composite moment-free Filon-type methods for the oscillatory integrals 
with smoothing integrand and without a stationary point. In sections 4 and 5, we develop the composite 
moment-free quadratures defined on a mesh according to the wave number k and the properties of the 
integrand for the oscillatory integrals with both singularities and stationary points. The formulas proposed 
in Section 4 have a polynomial order of convergence, and those in Section 5 have an exponential order 
of convergence. Numerical experiments are presented in Section 6 to confirm the theoretical estimates on 
the accuracy of the proposed formulas. Moreover, we compare the numerical performance of the proposed 
formulas with that of those recently proposed in mm- We summarize our conclusions in Section 7. 

2 The Filon-type Quadrature Method 

The goal of this paper is to develop quadrature methods for evaluating oscillatory integrals in the form 

Uf,g\:= dx, (2.1) 

where I := [0,1], k 2> 1 is the wave number, / € T 1 (J) has weak singularities and the oscillator g G C°°(I) 
has stationary points. Our main idea to fulfil this may be described as follows. We first develop a basic 
quadrature formula for computing an oscillatory integral defined on a subinterval [a, b] of I whose integrand 
has no singularities or stationary points. We then design an appropriate partition 0 = xq < x\ < ... < 
x n -\ < x n = 1 of I according to the wave number k, the singularities of / and the stationary points of g 
and employ the basic quadrature formula for each of the integrals defined on the subintervals [xj,Xj+ 1 ], 
for j € Z n _i := {0,1,..., n - 1}. 

We recall the Filon-type quadrature method proposed in m for computing the integral 

lt b] [/,</]== I'f(x)e iK9{x) dx, (2.2) 

J a 

where [a, b\ C I, f is continuous on [a, b] and the oscillator g is continuously differentiable on [a, b] and has 
no stationary point in [o, 6]. By a change of variables g{x) —>• x, the integral in (12.21) may be written as 

l [ : M lf,9}= [ 9(b) *(x)e iKX dx, (2.3) 

Jg(a) 

where 4/(a:) := ((// g') o g _1 ) (x), for x € [g(a),g(b)}. For a fixed positive integer m we approximate T by 
its Lagrange interpolation polynomial of degree m. Since g is differentiable on [a, b\ and has no stationary 
point in the interval, g must be strictly monotone on the interval. Without loss of generality, we assume 
that g is strictly increasing and otherwise we consider the integral [/, —g\. Letting m € N := {1, 2,...}, 
we choose m + 1 interpolation nodes in [a, b] such that a = to < t± < ... < t m -\ <t m = b. Let p m denote 
the Lagrange polynomial which interpolates T at the points g(tj), j € Z m . Hence, its Newton form is 
given by p m = a j w j with w j ( x ) = I1 igz,-i( x - 9(U)) for x € [g(a),g(b)\ and j € Z m , where the 

coefficients aj, j G Z m , are the divided differences of T, that is, aj := ^[g(to), g(t\),..., g(tj)\. When 
computing a,j, we are required to compute the values of T at g(tj). Noting that (g -1 o g) (tj) = tj for 
j G Z m , we need only to evaluate the functional evaluations of / and g' at the points tj. There is no need 


to calculate the inverse of g. A Filon-type quadrature rule is then obtained by replacing T in (12.31) with 
p m . That is, we use 

eL“l3[/,9]~ £ “A b(a) ' 9WI [“'j.9l. ( 2 - 4 ) 

j 

where g(x) = x for x G [g(a), g(b)], to approximate the integral (12.31) . In formula (|2.4j) . the integrals 
j-lFF-riF)} [wj,g\ for j G Z m can be computed exactly and efficiently. We postpone the description of 


computing these integrals and turn our attention to error analysis of formula 


For a function (j> G 


C(Q), let \\(/)\ 

^K,m[f,g] ■= ^’ b] [f,g] - QL“m [f,g\ and a := || g‘ 


:= max {|0(x)|}. Note that in this paper 12 will be either [a, b] or [g(a), g{b)]. Let 


estimate 


. According to m, we have the following error 

gfcg[/,d < 3(m + 2 1} T( m+1 ) a m (b-a) m . (2.5) 

’ m\n z ~~ 


Estimate (12.51) demonstrates that the decay of the error of the Filon-type method (12.41) is of order 0(k~ 2 ). 
The decay of (12.41) is also 0(k~ 2 ) for the linear oscillator (see, [HI .U, S3)- For a quadrature formula F 
that approximates integral (12.31) . we use N(F) to denote the number of evaluations of the integrand 4/ 
used in the formula. According to (12.41) . we have that M [o^n’mif, g]^j < m + 1. 

From (|2.5j) . we see that convergence of the quadrature rule Qjk’m[f,g] is affected by a. When a = 1, 
the error bound in (12.51) is not affected by a at all. When a < 1, the error bound in (12.51) decreases 
exponentially with respect to a as m —>• oo. When a > 1, the error bound in (12.51) grows exponentially 
with respect to a as m —>• oo. For the case a > 1 to accelerate convergence we subdivide [a, b] into N G N 
equal subintervals and approximate 4/ by its Lagrange interpolation polynomial p m of degree m G N on 
each of the subintervals. For IV G N, we use yj := a + (b — a)j/N for j G Zjy to denote the partition of 
[a, b] and [/, g\ computed by (|2.4I) to approximate [/, g\ for j € Z)y. This leads to the 

quadrature formula m [/, g\ for computing (12.31) , defined by 


QiiJhg] ■■= E 2fe 1,wJ [/.d 




( 2 . 6 ) 


Note that Qi'^ m [f,g\ = SL“m [/, g]. This quadrature formula will be used in the following sections for 
computing the oscillatory integrals with the integrand without a singularity or a stationary point. 

In the following theorem, we analyze the error £ [ ^ m [f,g] ■= T^' b \f,g] - • 


Theorem 2.1 If 'L € C m+1 [g(a), g(b)\, then for N G N 


c[a,b\ 

°N,K,m 


[f, g] < 


3 (m + 1) 
m!K 2 AC m_1 


Vj/lm+l) 


a m {b-a) r 


(2.7) 


and J\f ^Qjv’k ml/’^]) — N m + F I n particular, if N is chosen as \a], the smallest integer not less than 
a, then 


£ 


Ml 

N.K.m 


[f,g] < 


3(m + 1) 
min 2 

?[yj- iiS/j] 


a (p _ a )"*. (2.8) 

OO 

Proof: Let N G N. We estimate the error 1,WJ [f, g] ■= [/, g] — Qfff m 1,v ^ [/, g] by first em¬ 

ploying estimate (|2.5I) with [a, 6] replaced by [((/j_i, yf and then Summing up both sides of the resulting 
inequality over j G Noting y-j — yj-i = (b — a)/N , this leads to estimate (12.71) . According to the 

algorithm (12.61) . noting that the nodes yj for j G Z)^_ 1 are used twice in the algorithm, we conclude that 

J\f (QJv’k m [/’ d ]) — -AWA (Si"fm/iS]) — (N — 1) < Nm + 1. In particular, when N is chosen as [cr], we 
substitute N = [a] into estimate (I2.7j) to yield estimate (12.81) . □ 

























Comparing estimate (12.71) in Theorem 12.11 with estimate (12.51) . we see that formula (12.61) uses Nm + 1 
number of the functional evaluations of T, which is N times of that used in the traditional Filon method to 
reach the order of error estimate N~ m+1 £^m [/, g\- Formula (12.61) will serve as a basic quadrature formula 
in this paper for developing sophisticated formulas for computing singular oscillatory integrals. 

In the remaining sections of this paper, we shall consider the following three cases: 

(i) When / and g are smooth and g has no stationary point or inflection point in I , according to 
the wave number k we design a partition 0 = xq < x\ < ... < x n = 1 , and write I K [f,g\ = 

[f,g]- Formula (12.61) is then used to compute integrals Z^ 3 ~ 1,Xj \f, g\ for j G Z+. 

(ii) When / has a weak singularity only at the origin and g is smooth without a stationary point or 
an inflection point in I, we first divide / into two subintervals [0, 6 ] and A := [6,1] such that the 
integrand of zl°’^[/, g] does not rapidly oscillate and that of l£[f,g] does not have singularity. The 
integral zJ-°’^ [/, g] is calculated by a quadrature rule for weak singular integrals. The integral Z^[/, g] 
is computed by the method described in item (i). 

(iii) When / has a weak singularity only at the origin and g is smooth with one stationary point at the 
origin and has no inflection point in I, we first divide I into two subintervals [0, 6 ] and A such that 
the integrand of Z^ 0,6 ^ [/, g] does not rapidly oscillate and that of Z^[/, g] does not have singularity or 
stationary point. The integrals zJ-°’^[/, g] and Z^[/, g] are handled in the same way as (ii). 

The case of / having a finite number of singularities in /, g having a finite number of stationary points 
or inflection points in / can be treated by splitting / into subintervals, on each of which / has only one 
singular point or g has only one stationary point at an end-point and without an inflection point. 

3 The Composite Filon-type Quadrature Method 

In this section, we develop composite Filon-type qadrature methods for computing the oscillatory integrals 
m- Specifically, we assume that f,g£ C 2 (I ) and g is increasing monotonically with ||g'|| = er, 

g'{x) 7 ^ 0 for x £ I and g"{x) / 0 for x G (0,1). We shall partition the interval I into n subintervals 
according to the wave number k, and propose a composite quadrature rule, where we approximate T by its 
Lagrange interpolation polynomial of variable degrees in each subinterval, aiming at the asymptotic error 
order (D{K~ n ~ 1 ). 

We first motivate the construction of a K-dependent partition of I. By a change of variables of kx —>■ x, 
the integral ( 12 . 11 ) becomes 

I K [f,g] = I/* £ f(x/ K )e iK ^dx + l/ K £ f(x/K)e iK ^dx. (3.1) 

The integrals on the right hand side of (13.11) do not oscillate rapidly since | (ng(x/ k))'\ < a for x G [0, k\. 
However, traditional quadratures for computing the second integral on the right hand side of (13.11) lead 
to prohibitive costs for a large k. Inspired by the quadratures for singular integrals using graded points 
proposed in m, for n G N with n > 1 we suggest the partition of [1, At] with graded points l ^/( n 
for j G Z+. Using the change of variables x/k —>• x, we obtain a partition of I. 

We now describe the construction of the K-dependent partition of I. To this end, we fix k. For n G N 
with n > 1 , let n K denote the partition of / with nodes defined by 

xq = 0, Xj = k^ _ 1 )/(” _1 ) _1 , for j G Z+. (3.2) 

According to the partition n K , the integral (12.11) is written as I K [f,g] = yT eZ + zjf^ -1 ’ 3 ^ [/, g\. Computing 

the integral Z K [f,g] is then reduced to computing the integrals [/, g] for j G Z+. We shall use 

formula (12.61) to calculate these integrals. For this purpose, we define the quantities 

Mj := max{|^'(xj_i)| , \g\xj)\} and Nj := \M j ], for j G Z+. 


Since g is monotonically increasing on /, we have for j £ Z+ that 

Mj = max{|</(x)| : x £ [xj_i,Xj]} and Mj < a. (3-3) 

We shall develop two quadrature methods. Method one uses a fixed number of quadrature points in each 
of the subintervals and has a polynomial order (in terms of the wave number) of convergence. Method two 
uses variable number of quadrature points in the subintervals and achieves an exponential order (in terms 
of the wave number) of convergence. 

We first describe the method having a polynomial convergence order. We choose a fixed positive integer 
m. For each j £ Z+, we use [/, g\ defined as in (|2.fil) to approximate [J, gf Integral 1 K [/, g] 

defined by (12.ID is then approximated by the quadrature formula 

Qre,n,m [/j S'] : = ^ [/) d\ • (3-4) 

j 6Z+ 


In the next proposition, we provide an estimate of the error £ K , n ,m.[f, g\ ■= |T h ,[/,g] ~ QK,n,m[f, g]\ and the 
number Af (Q Kt n,m[fi d]) of functional evaluations used in the formula Q K} n,m[fi9\- 


Proposition 3.1 For with n > 1. let g := max 

some mgN, then 


£K,n,m[f j g\ F 


3 (m + 1) 
m\n 2 


{1/k,1 

^(m+l) 


i)}. p e C m+1 [g(0),g(l)} for 


arj' 


jm— 1 


and, Af (Q«, n , m [/, g]) < \a]nm + 1 . 


Proof: The proof is done by applying Theorem 12.11 on each of the subintervals [xj-i,Xj\. Specifically, 


for j £ Z+, we use Theorem O to estimate £j [/, g] := 


i [ : i - 1 ’ Xi] [f,g\ - Q [ N j r^ ] [f,g\ 


For j £ Z+, 


we apply (12.81) with a being replaced by Mj and b — a by hj := Xj — Xj -1 to conclude that £j [/, < 7 ] < 
||^r( m+1 )|| oo Mjh™. Note that h\ = xi — xq = 1/k < g, and for j € Z+ with j > 1, hj = 


Xj — Xj -1 = Xj (l — < Xjr/ < 77, and Mj < < 7 . Substitute these bounds into the inequality above 

and summing the resulting inequality over j £ Z+, we obtain the desired estimate for £ K ,n,m[f, g\- 

Using Theorem [2J] again yields that Af (Q K ,n,m[/, g]) < J2j^z+ (Afjm + 1) — (n — 1) < [~cr]nm + l. □ 


As a direct consequence of Proposition 13.11 we have the next estimates for the case having the linear 
oscillator g(x) = x, x £ /, where a = 1 . 


Corollary 3.2 For n £ N with n > 1, let g := max {1/k, 1 — k 1 ^ n T}. If f £ C m+1 (I) for some m £ N 
and g{x) = x for x £ I, then £ K , n ,m[f,g\ < ||/ (m+1) IL V^ 1 and Af (Q K ,n,m[/, g}) <nm + 1. 

TTl'.K, 


We now turn to developing the quadrature formula having an exponential convergence order. This is 
done by choosing variable numbers of quadrature nodes in the subintervals of I. Specifically, for n £ N 
with n > 1 and for the partition of I chosen as (13.21) . we let mj := \n(n — 1 )/(n + 1 — j)~\ for j £ Z+. For 
each j £ Z+, we use [/, g] to approximate l]f j ~ 1,x ^[f,g\. Integral Z K [/, g] defined by (12.11) is then 

approximated by the quadrature formula 

0n,n[f,9]-= E (3-5) 

jezt 


We next study the error £ K)U [f,g\ := \I K {f,g\ - Q K , n [f,g}\ of the quadrature Q K , n [f,g\ ■ To this end, we 
first establish two technical lemmas. 


Lemma 3.3 There exists a positive constant c such that for all n £ N with n > 2, z + 1 /{mj — 1) < c. 












Proof: We prove this result by estimating the lower bound of the set {rrij — 1 : j G Z+}. For j G Z+, we 
have that rrij — 1 > n(n — l)/(n + 1 — j) — 1 > n — 2. Thus, we obtain that V( m j — 1) < n/(n — 2 ), 

which is bounded by a constant c. □ 


Lemma 3.4 There exists a positive constant c such that for all k > 1, n G N that satisfy 


(n — l)(ln (n — 2 + e) — 1) > InK, (3.6) 

and j G Z+ with j > l, 

( K V(™-i) _ i) m ^- 2 /(m i - 2)! < c(n - 2)~ 1/2 . 

Proof: By the definition of mj, we see that rrij — 2 > n — 2 for j > 1. Condition (13.611 implies that n > 2 
since k > 1. We observe from the Stirling formula [T] for n G N 

n! > y/2rm (n/e) 77 . (3.7) 

Using inequality (13.711 with n := TOj — 2, we conclude that there exists a positive constant c such that 
for all n G N with n > 2 and j € Z+ with j > 1, \/{mj — 2)! < c{mj — 2 ) -1 / 2 (e/(m,j — 2)) mj ~ 2 < 
c(n — 2 ) -1 / 2 (e/(n — 2)) mj ~ 2 . On the other hand, condition (13.611 implies that e^ 1 ^ 71-1 ) — 1) < n — 
2. This together with the above inequality ensures that there exists a positive constant c such that 
for all k > 1, n G N satisfying (|3.6I1 and j G Z+ with j > 1, (^VO -1 ) _ \) m i~ 2 /{mj — 2)\ < c(n — 
2 ) -1 / 2 (e(«: 1 /( n_1 ) — 1 )/(n — 2)) mj 2 < c(n — 2)” 1 / 2 . This concludes the desired result. □ 

For a function f G C' oo (0), we let ||</>|| n := max {||(/>^)|| : j G Z n } for n € N. We are now ready to 

establish the estimate for £ K>n [f,g\. 


Theorem 3.5 If 4' G C' oo [< 7 ( 0 ), g(l)], i/ien there exists a positive constant c such that for all 1 1 > 1 and 
n G N satisfying m, 

£*,»[/> 9} < ca{n - 2)~ 1 / 2 K~ n ~ l ||T || (n _ 1)n+1 . 

For n G N with n > 2, there holds the estimate Af (Q K ,n[fi 5]) < f er~| (n(n — 1) In n + n 2 ) + 1. 


Proof: We establish the error bound by estimating errors £j[f,g] ■= T/f 3 g\ — Q}/l [/, g\ for 

j G Z+, and then sum them over j. By employing (12.811 with b — a being replaced by hj and m by mj and 
(S3]), we obtain for j G Z^ that 




mj\n 2 ~ 3 




(3.8) 


For j = 1, we have that £\ [f,g] < 
into (|3.8|) we obtain that 


6 an 


—n —1 


|'p( n )|| . For j > 1, by substituting hj = k **-1 (ac 1 ^ 71 — 1) 


(77.-2)! 


£j[f,g\ < 


6 a ( k 1 /^-!) _ 1)^-2 

mj — 1 (mj — 2 )! 


,l/(n-l) _ -O m i 2 ^(rrtj+l) 


(3.9) 


On the other hand, condition (13.61) implies that n > 2. Thus, (ac 1 ^ 71 ^ — l)“ < k. Applying Lemma [37 
to (|3.9p yields a positive constant c such that for all /c 1, ti G N satisfying m, and j G Z+ with 

j > 1, £j[f,g\ < mj ~ 1 || v I' (mi+1) || 00 - Since for j € Z+ with j > 1, J n 


m 9 - — 1 < 
71-1 3 


j — n — 1 n{n — 1 ) 


n — 1 n + 1 — j 


7 — 1 = —n — 1 , we conclude that there exists a positive constant c such that for all k > 1 , 


n G N satisfying (13.61) . and j G Z+ with j > 1, [/, g] < ca 


- 2 W 1 / 2 


in- 2 ) 


- 1 ||\J/( m j+ 1 )|| . Summing 


up 


mj - 1 























the bound of errors Sj [/, g] over j € Z+ and applying Lemma 13.31 we observe that there exists a positive 
constant c such that for all k > 1 and n € N satisfying m 


6a 


£ K Af,9] < 


—n— 1 


_n_r I A(n-2)-V 2 

'L (n) + C<7K n 1 < > ---- 

™ l IB; — 1 


I'Ll 


(n— l)n+l 


. 1=2 


< ca(n — 2) V 2 « n 1 |l^ll(n-i)n+i • 


It remains to estimate the number of functional evaluations used in the quadrature formula. To this 
end, we note that Theorem 12. II yields 

A/" (Q K ,n[f,g}) < { T 0- ! (n(n-l)/(n + 1 - j) + 1 ) + i} - (n - 1.) 

< [cr] < n + n(n - 1 ) y] l/(w + 1 — j) > + 1 . 

I i£Z+ J 

For n € N by using X^<=z+ 1/j < hrn + 1 we have that 

A f(Q K , n [f,g]) < H {n + n(n - l)(ln(n) + 1)} + 1 = \a] (n(n - l)lnn + n 2 ) + 1 , 
which completes the proof. □ 


As a direct consequence of Theorem 13.51 we have the next estimates for the case having the linear 
oscillator. 


Corollary 3.6 If f € then there exists a positive constant c such that for all k > 1 and n € N 

satisfying (|3.6I) . £ K ,n[f,g\ < c(n — 2)~ 1 ^ 2 n~ n ~ l ||/||( n _ 1 ) n+1 - For n € N with n > 2, there holds the estimate 
A r(Q K ,n[f,g]) < n(n- 1 ) In n + n 2 + 1 . 

The Filon-type methods mm and the Levin-type methods EU|29] achieving a convergence order 
higher than 0 (k _ 2 ) require the evaluation of derivatives of / and g'. The Filon-Clenshaw-Curtis rules jT] 
for the oscillatory integrals with nonlinear oscillator requires the evaluation of g~ l . Quadrature methods 
developed in this section do not require computing derivatives of / or those of g' , nor evaluating g v . 
Moreover, since condition (13.61) implies that InK < nlnn, Theorem 13.51 demonstrates that the quadra¬ 
ture Qn,n[fig\ achieves the asymptotic convergence order 0 (n~ 1 // 2 K _Tl ~ 1 ) with only 0( In 2 k) number of 
functional evaluations. 


4 Quadratures with a Polynomial Order of Convergence 

In this section, we consider computing the oscillatory integral (ED), where / is allowed to have weak 
singularities, g has stationary points and has no inflection point in I. We shall develop quadrature formulas 
having a polynomial order of convergence. 

The key idea to be employed is to split the interval / into two subintervals such that on one subinterval 
the integrand has a weak singularity but no oscillation and on the other subinterval it has oscillation but no 
singularity or stationary point. We then employ quadratures using graded points for the singular integral 
and design a composite quadrature rules using a partition of the subinterval which is formed according to 
the wave number k and the property of g for the oscillatory integral. 

We first describe the weak singularity of a function defined on I according to m- Let 5 be a subset 
of I containing a finite number of points. For some a G (—1,1) and a nonnegative integer m, a real-valued 
function / € C m (I\S ) is said to be of Type(a, m, S) if there exists a positive constant c such that for all 
x € I\S, 


f im) (x) 


<c[us(x)] a - m 


(4.1) 










where the function us associated with S is defined by us(x ) = inf {\x — t| : t G S} for x € I. When we 
say that / is of Type(a, oo, S), we mean that for all m € N, / S C°°(I\S ) satisfies (14.11) . The parameter 
a is called the index of singularity. For a > 0, this notation was introduced by Rice in [32]. 

In this section, / is allowed to have a single weakly singular point at zero with index p, € (—1,1) and 
g satisfies the following assumption: 

Assumption 4.1 For a nonnegative integer r, the function g € C r+1 (I) has a single stationary point at 
zero satisfying g^\0 ) = 0 for j € 7L r , g^ r+l \x) ^ 0 for x € I, and a(r) := ||g( r+1 ) 11/(7" + !)■ ^ an d 9 
is increasing monotonically without an inflection point in I. 

The requirement that g(0) = 0 in Assumption 14.11 is without loss of generality, since if g(0) ^ 0, we 
compute instead the integral Z K [/, g — g(0)] + exp {iKg(0)}l K [f, 0]. For the case r = 0, the oscillator g does 
not have a stationary point. 

We now write the integral p.ll) as the sum of two integrals: a weakly singular integral without rapid 
oscillation and an oscillatory integral without a singularity or a stationary point. According to the as¬ 
sumption on g, by the Taylor theorem, for each x € / there exists a constant f x € [0, (Kcr(r))~ 1 ^ r+1 'lx\ 
such that k\ g [(Ka(r))~ 1 ^ r+1 ' > x) | = \g tyr+1 \f x )x r+1 / a(r)\/ (r + 1)!. Hence, we see that n\g(x)\ < 1 for 
0 < x < (n(j(r))~ 1 K r+l \ Thus, for such an x the function e 1K does not oscillate rapidly. Motivated 
from the above discussion, we introduce 


^cr(r) 


K, cr(r) < 1 , 

k a(r), a(r) > 1 , 


and define 


X r := K 


-l/(r+l) 

(r) 


(4.2) 


We split the interval I into two subintervals [0,A r ] and [A r ,l]. Correspondingly, integral ( 12 . 11 ) may be 
written as the sum of integrals on these two subintervals. Let A := [A r ,l] and for (f> € L l (I) we set 
m := fo (f(x)dx. Using a change of variables: y = A^ 1 ^ for x € [0, A r ], the integral (| 2 . 1 I) is rewritten as 
Zk [f, g] = 2 T [/, 9 ] + [/> 9\ ) where 

l»[f,g]:=\rI[<p K ], (4.3) 


with 


tp K (x) := f (X r x)e iK ^ XrX \ for x G / 


(4.4) 


and 

l£[f,g]:= £ f(x)j*3Wdx. (4.5) 

Note that the function tp K defined as in (I4.4|) has a singularity at the origin but has no oscillation and 
the integrand in (|4.5|> has no singularity or stationary point but has oscillation. We shall treat these two 
integrals separately. 

We shall develop a quadrature rule for computing the singular integral (14.31) having a polynomial order 
(in terms of the number of nodes used in the partition) of convergence, and a quadrature rule for computing 
the oscillatory integral (14.51) having a polynomial order (in terms of the wave number) of convergence. 

We now consider the integral I[(p K ] that appears in (14.31) . The integrand (p K defined by (14.41) does not 
oscillate rapidly. The classical quadrature rules for weakly singular integrals developed in m can then be 
used to treat the singularity. Below, we briefly review the quadrature rules. 

We begin with describing the Gauss-Legendre quadrature rule for integral [tjj] := f b if>(x)dx, where 
if is a smooth function defined on [ a,b}. Given m € N, we denote by — 1 < t\ < t 2 < ■ ■ ■ < t m -i < 
t m < 1 the zeros of the Legendre polynomial P m of degree m and by Uj := 2(1 — tj)[mP m -i(tj)}^ 2 

for j € Z+, the weights of the Gauss-Legendre quadrature rule which has the form [if] := ( 6 - 
a )/2 ([(& — a )tj + (b + a)\/ 2). There is a constant £ e [a, b\ such that the error of the approxi¬ 
mation to Z [“ ,6 1 [if] is given by lZm’ b \if] := (b — a)‘ 2m+1 if^ 2rn \f) /( 2 2 m ( 2 m + 1 )!). 


We now recall the integral method proposed in m ■ Given m G N, let p := (2 m + 1)/(1 + /i). For 
s G N with s > 1, according to the parameter p we choose s + 1 points given by Xj := s~ p j p for j G Z s . 
The quadrature rule for 1[ip K ] is obtained by replacing ip K on [xq,x\] by zero and using Qm J ’ Xj+1 ^ [p K ] for 
computing the integrals I^ x ^ x i+^[ip K \ for j G Zj~_i- Integral X p [f,g\ defined by (I I.-I I) is then approximated 
by the quadrature formula 


■= X r 


E 


a 


X jl X j + 1 ] 


Wk 


(4.6) 


We next estimate the error £^ m [f,g] '■= | — 2T [/, g] |. We need a lemma that estimates the 

norm of w K (x ) := exp{u K (x)}, for x € /, where u K (x ) := i Kg ( \ r x ) for x G /. This requires the use of the 
Faa di Bruno formula [IS, 30] [33] for derivatives of the composition of two functions. For a fixed n G N, if 
the derivatives of order n of two functions f> and if are defined, then 


(4>o -0) (n) = ^2 • • • ,^ (n J+1) ) , 

jezi 


(4.7) 


where for j G Z^, 


Hn,j (•Z'l) X 2; ■ ■ ■ ; j • 1 ) — ^ ) 


n! 


mi\m 2 \ ■ ■ ■ m n - j+ i\ 


n * 




Xl\ m l 


l\J 


(4.8) 


where the sum is taken over all (n—j+l)-tuples (mi, ..., m n _, + i) satisfying the constraints ~+ m l = 
j and ^2i £Z + Ifni = n - For n G N and for j G Z+, we let B n j := B n j( 1,1,... , 1) for brevity. Note that 
B n j = 1/j! 'Y^i &z X—X)i- l C l jl n are the Stirling numbers of the second kind [1]. They have the property m 
that for n € N with n >2 and j G Z+ satisfying j < 2 _1 n, B n ^ n _j + \ < B n j. For n G N, the Bell number 

/ q \ ^ 

has the bound [3j that B n ^2j eZn B n j < ^ + 1)/ ' 

Lemma 4.2 If g G C 2m (I ) /or some m G N satisfying Assumption 14.11 f/ien f/iere exists a positive constant 
c such that for all k > 1, ||u> K || 2m < c. 

Proof: For n G Z^" m , by using (14.71) . we obtain that 


«>l n) = J+1) ) • 

jez+ 

From the definition of u K , we have two inequalities below. If j G Z+ with j > r + 1, we have that 


(4.9) 


u^\x) 


< n\ J r 


g W (Ar*) 


< 


for x G /. 


For j G Z+ with / < r + 1, by Assumption 14.11 there exists a constant f x G [0, A r x] such that 

g( j) (X r x) = g^Xix)/{r + 1 - j)! (A r x) r+1 ^ , for x G /, 
which implies that for all j G Z+ with j < r + 1 and iGl, 


(4.10) 




< kA^ 


{Kx) 


— Ilffllr+l ' 


(4.11) 


Letting M := max {1, ||fif|| r+1 
and x G 7 that 


12m 


n G Z+ m that 


This together with 


B n j ( ,. 


(n) 


} and applying (14.811 with (14.1011 and (14.111) . we observe for all k > 1 
< B n jMB This together with (14.911 yields for all k > 1 and 


(n-j+l) 


OO 

W, 


< KILE,-^ Bn.jA'B < II^kIIqo B n M n , which is a constant independent of k. 


^ II OO 


< 1 yields the desired estimate. 


□ 





















We next estimate the derivatives of the function tp K defined by (I4.4[) . We recall the Leibniz formula for 
the n-th derivative of the product of two functions for n € N, 

OV0 (n) = E (4.12) 

mEZ n 

Tl\ 

where C™ ——-— are the binomial coefficients that satisfy 

m!(n — m)\ 

E C n= 2 n - (4.13) 

mez„ 


Lemma 4.3 Let m € N. If f is of Type(/r, 2m, {0}) for some p € (—1,1) and g 
Assumption 14.11 then there exists a positive constant c such that for all re > 1 and x € 


€ C 2m (I ) satisfies 


( 0 , 1 ], 


(2m) 


(x) 


< 


C \rX tl ~ 2m . 


Proof: Applying the Leibniz formula (I4.12|) to the function yields 


4 2 m) W= E CL(/(A r x)) U) w< 2 ”" fl W, for * 6 ( 0 , 1 ]. 


(4.14) 


From the assumption on /, there exists a positive constant c such that for all re > 1, j G Z 2 m and x € (0,1], 
(/(A r x))^ < |/0) (A r x)| < cXrX^-i. This together with ()4.14j) yields a positive constant c such that 

'or all re > 1 and x € ( 0 , 1 ], 


4 2m, w| < c E cLk^-’ 


w 


(2m-j) 


(x) 


< cA(fx M_2m 


re 


K 112m ^ . ('2;n • 

j 


Using Lemma 14.21 and formula ()4.13j) in the inequality above, we obtain the desired estimate. □ 

We need a technical result for the integral of a function of Type(^, 0, {0}) for some p G (—1, 1). 

Lemma 4.4 If <f> is of Type(/i, 0, {0}) for some p € (—1,1), then there exists a positive constant c such 
that for all 0 < b < 1 and g > 0 , Jq \4>(gx)\ dx < cg /1 b 1+/J ’. 

Proof: We prove this result by bounding \4>(gx)\ by cg^x^ 1 and computing the resulting integral exactly. 

□ 


For a quadrature formula F for approximation of (|4.3I) . we use Af(F) to denote the number of evalua¬ 
tions of the integrand (p K used in the formula. With the above preparation, we estimate the error £^ m [f, g] 
and Af (Q^ m [/> g\) i n the following theorem. 

Theorem 4.5 Let m € N. If f is of Type(p,2m,{0}) for some p € (—1,1) and g G C 2m (I ) satisfies 
Assumption 14.1L then there exists a positive constant c such that for all re > 1, s G N with s > 1 


9 ] <« 


(4.15) 


and AT {Q S fj, m [f, g]) < (s — 1 )m. If g{x) = x for x G I, then the upper bound in (14.151) reduces to 
cre" 1 "^" 2 ™. 


Proof: Let £ 0 [<£>«) := <p K (x)dx and £j[<p K \ := Qm ’ Xj+l1 YPk] ~ X^' x i+^ [<p K ] for j G Z+_ 1 . We prove 
(14.151) by estimating Sj\<p K \ and then sum them over j G Z s _i. 

We first consider £o[</j k ]. By applying Lemma [4.41 with g := X r and b := x\ = s 
conclude that there exists a positive constant c such that for all re > 1 and s G N with s > 1 


X 1 (2m+l)/(l+/d 


we 


€0 



—2m—1 


|/ (A r x)| dx < cA(ls 


























Note that X r is defined by (14. 2\ . For j G Z+_j, by the error bound of the Gauss-Legendre quadrature, there 


exists aconstant € [xj,Xj+ 1 ] such that £j[<p K \ < h^ m+L /( 2 2 m ( 2 rn+l)!), where hj := Xj + \—Xj. 

Using Lemma 14.31 with fj > xj for j G ZjLi, there exists a positive constant c such that for all k > 1, 

Substituting this inequality into the inequality 
obtained from the proof of Theorem 2.3 in |j l9l . 


< c\rXj 


—2m 


—2m— 1 


s € N with s > 1 and j G m (£j) 

above and employing the estimate h? m+1 Xj 2m < cs 
we see that there exists a positive constant c such that for all k > 1, s € N with s > 1 and j € Z+_ l5 
< cXr s~ 2m—1 . Summing up the bound of errors over j € Z s _i and using the definition (14.21) 

of X r we obtain the estimate (|4.15|) . 

The bound on the number of functional evaluations used in the algorithm (|4.6I) may be obtained directly. 
When g(x) = x, we substitute r = 0 and k ct (o) = n into estimate (|4.15l) to yield the special result. □ 

We next develop a quadrature formula for computing the oscillatory integral (14.51) . For n € N, we 
choose the partition II K ^( r ) of A with nodes defined by 


T . — K U/n-l)/(r+l) ■ j 

x 3 •— K cr(r) > 4 t ^ n - 


(4.16) 


According to the partition, the integral (14.5p is rewritten as X£[f,g] = ljf j 1 ,Xj \f,g\. For each 

j € Z+, we use quadrature formula (12.61) to compute an approximation [/, g] of the integral 

[y,^], w ]- iere _/v) is a positive integer to be specified later. We estimate the error £j[f,g\ := 
jftj-uxj .] [/, g] _ g\ for j € Z+ in the following lemma. To this end, we define 0(y) : = 

^(g(l)y) for y G /, d(r) := min {|</ r+ 1 )(x)| /(r + l)!:iG /}, /3 n (r) := X — 1 with A r defined by (I4.2p . 
Mj '■= max {\g\xj-i)\ , |fl / (x j )|} and qj := MjXj-i/g{xj-\). for j G Z+. 

Lemma 4.6 If 0 is of Typ e(a,m + 1, {0}) for some a G (—1,1) and some m G N, and g satisfies 
Assumvtion 14.11 then there exists a positive constant c such that for all k > 1. n G N and j G Z+ 


£j[f,g\ < c 


isoor 


(Mr))’ 


(m — 1 )\k 2 N, 


Proof: We prove this result by applying (12.71) in Theorem 12.11 to £j[f,g] for j G Z+. From (12.71) . we have 


for all j G Z+ that £Af,g\ < 3 ^ m + 1 ^-, ||T( m+1 )|| MW, where hj : 

J n juizn — o AT-m-1 II Moo 3 3 ’ ■> 


m\n 2 N, 


= x 7 — Xj_i. In the estimate 


above, using the relation vfd m+ 1 )(< 7 (l)?/) = (< 7 ( 1 )) m !0( m +i)(y) for y G (0,1] with the assumption on 0, 
we observe that there exists a positive constant c such that for all k > 1, n G N and j G Z+, 


[./) 9] < c 


b(l)l' 


(m — l)\n 2 N r - 


-q™(g( Xj . 1 )) 


a-l -ro 7 m 


Note that (g(xj_i)) Q 1 < (<5(r))“ 1 x ( j'f h 1 1 ^ a ^ and Xjf n 1 hJ L = (/3 n (r)) m for j G Z+. Substituting these 
relations into the inequality above yields the desired result. □ 


We estimate an upper bound of the quantity := max )gZ + {dj}- 

Lemma 4.7 If g satisfies Assumption ^ TJ then for all n G N and k > 1, Q n < (r + l)K r J^'a(r)/5(r). 

Proof: Using the definition of 6(r ) and <r(r), we have that g(x) > 5{r)x r+l and |< 7 / (a:)| < (r + l)<r(r)x r for 
x G I. Thus, \g(xj-{)\ > 8{r)x r £f\ and Mj < (r + l)cr{r)xl. For n G N, we obtain for k > 1 and j G Z+ 

□ 


that qj < (r + l)^" (r+ 1 ) V(r)/d(r). 















We now discuss the choice of Nj. Lemma 14.71 demonstrates that the upper bound of Q n is independent 
of k for the case r = 0 and it depends on k for the case r > 0. Therefore, we need to consider these two 
cases separately. For the case r = 0, we choose 


Nj ■■= I" Qj] for j € Z+. 


(4.17) 


For the case r > 0, we choose 


Nj-.= 


ml (m— 1) 
3 


for j € Z,1 


(4.18) 


so that Nj (rn ' qj 1 ' for all j € Z+ are independent of k. With this choice of Nj, we use ^\f,g] to 

approximate the integral zjf[/, g] for j € Z+. Thus, integral T^[f, g\ defined by (14.5|) is approximated 
by the quadrature formula 

Qin,m[f,9] ■= E (4-19) 

jGZin 

We next estimate the error £^ n>m [f,g] := \ 1^[f,g\ - Q^ n>m [f,g]\ and Af g}) ■ We first con¬ 

sider the case r = 0. For simplicity of notation, we let 5q : = 5(0), Uq := <x(0) and := j3 n { 0). 


Theorem 4.8 Let r = 0. If 0 is of Type(o, m + 1, {0}) for some a € (—1,1) and some m G N and g 
satisfies Assumption 14.11 with r = 0 , then there exists a positive constant c such that for all k > 1 and 
n € N satisfying n > In K ao j In 2, 

£ K,n,m[f'9\ - c 5^- 2 a 0 K- 2 nlf a lnK ao (/3 °) m_1 . (4.20) 

There holds the estimate N {Q^ n m \f , g]) < \ao/5o]nm + 1. 

Proof: We prove (14.201) by estimating £j[f,g] for j € Z+, and then summing them over j. By employing 
Lemmas 14.61 and l4~7l with r = 0 and the choice (14.171) of Nj, there exists a positive constant c such that 
for all k > 1, n € N and j € Z+, £j[f,g] < c 15(1)1 Sq ^ 2 o'ok _ 2 x“T 1 1 Since — 1 < a < 1, we see 

that x°jZ\ < Xq 1 = K af a - Substituting this inequality into the inequality above yields that there exists a 
positive constant c such that for all k > 1, n € N and j € Z+, £j [/, g] < c\g(l)\~ a 5 q~ 2 n\~ a (/3^) m . 
Moreover, condition n > In n ao /In 2 implies that n/3 ° < 21nK (T0 . Thus, there exists a positive constant c 
such that for all k > 1 and n € N satisfying n > In k CT o / In 2, 

£ K,n,m if > ff] < E < c|ff(l)r Q 5o” 2cJ OK“ 2 ^““lnK (TO , 

which is the desired estimate (14.201) . 

The estimate of A f ( Q^.n m [/> d \) is obtained by using formula (j4T9j), the choice of Nj and Lemma l4~7l 
with r = 0 . □ 


In passing, we comment on the hypothesis imposed to 0 in the last theorem. If / £ C m (I) for some 
m € N and g € C m (I ) satisfies Assumption 14.11 then 0 is of Type(— r/(r + l),m, {0}). A proof of this 
conclusion may be found in [7j. As a direct consequence of Theorem 14.81 we have the next estimates for 
the case with the linear oscillator g(x) = x for x € /, where we have that g{ 1 ) = 5o = 0 o = 1 and K ao = n. 

Corollary 4.9 If f is of Type(a, m + 1, {0}) for some a G (—1,1) and some m € N, then there exists a 
positive constant c such that for all n > 1 and n € N satisfying n > In n/ In 2, £ff n m [f, g] < ck - " -1 In k#™” 1 , 
where 9 n := — 1. For n € N, there holds the estimate N < nm + 1. 


We now consider the case r > 0. 












Theorem 4.10 Let r > 0. If 0 is o/Type(a, m + 1, {0}) for some a G (—1,1) and some m G N and g 
satisfies Assumvtion \AA[ then there exists a positive constant c such that for all k > 1 and n G N satisfying 
n > \n.K a i r )/\n2, 

£ K,n,m[f’9\ < c(r + 1) \g(l)\~ a (<y(r)) a- 1 K - 2 /c^ r “lnK CT(r) (/3 n (r)) m_1 . (4.21) 


There holds Af (Q% m [f,g]) < |"(2(r + l)a{r)/8(r)) m/{m X) 


nm + 1 . 


Proof: Estimate (I4.21j) is proved in a similar way to that of estimate (I4.20|) in Theorem 14.81 The only 
difference is that instead of using the hypotheses and results for the special case r = 0 , we use those for 
the case r > 0. We leave the details to the interested reader. 

It remains to estimate M (Q« n m [/,<?]) • According to formula (|4.19|) and Lemma l4~7l we obtain that 


N [/;#]) < E jeZ + { N j m + !}-(«-!)< ((r + 1) <r(r)/S(r)^ 


mj (m— 1) 


nm + 1. Note that 


condition n > In K a ^ r ) / In 2 for n G N implies that 


1 In 
K , , 
cr(r) 


< 2. 


This concludes the last result. 


□ 


To close this section, we point it out that the oscillatory integral (|2.1|) may be computed according to 
the formula Q s K ff,m [/, g] := Q^mifid] + with two fixed positive integers rh and m. 


5 Quadratures with an Exponential Order of Convergence 

In this section, we develop a quadrature method for computing the oscillatory integral (12.11) having an 
exponential order of convergence. As in section 4, we shall write integral (12.11) as the sum of a weakly 
singular integral (14.31) and an oscillatory integral (14.51) . and treat them separately. Specifically, we shall 
develop a quadrature rule for computing (|4.3I) having an exponential order in terms of a constant 7 € (0,1) 
of convergence and a quadrature rule for computing (|4.5I) having an exponential order in terms of the wave 
number of convergence. 

We now describe the quadrature rule for the integral X[p K ] that appears in (14.31) . In developing the 
quadrature rule for computing Z[<£> K ] we borrow an idea from [35] (see also [4]). We begin with describing 
a partition of I. For 7 € (0,1) and s € N with s > 1, let II 7 be a partition of / with nodes defined 
by xq = 0 and Xj = 7 s- - 7 for j G Z+. For e > 1, we choose rrij := \je] for j G Z^"_ 1 . The quadrature 

rule for computing T[<p K \ is formed by placing tp K on [xo,a:i] by 0 and using Qm’ ,Xj+ 1 ^ [<p K ] for computing 
the integrals L^ Xj,Xj +^[(p K ], for j G ,. Integral X^[f,g] defined by (14.31) is then approximated by the 
quadrature formula 

Q s 1 [f,g\:=Xr £ g£f J+1 W (5.1) 

j&i-i 

where A r is defined by (|4.2[) . We next estimate the error £^[f, g] := |X ft [/, g] — Q^[f, g]\. To this end, we 
impose the following hypothesis on g. 


Assumption 5.1 There exists a positive constant ( > 1 such that for all n G N, ||g|| n < C- 

We first estimate the norm of w K (x) = exp{tt K (x)} with u K (x) = ing (X r x) for x G / as Section 4. 

Lemma 5.2 If g G C°°{I) satisfies Assumption [54) then for all k > 1 and n G N, ||u; K || n < B n (f 1 . 

Proof: The proof of this lemma is similar to that of Lemma 14.21 Using Assumption 15.11 and applying 
(14.8|) with (14.101) and (14.111) we observe that for all n G N, k > 1, j G Z+, l G Z + and x G I, 


B j,i 


(j'-i+i) 


< But 1 . 



















This together with wiP = Yhi^z+ w ^j,i ■ ■ ■ -,Uk l+1 ^ yields that for all n E N, j E Z+ and 


k > 1 


w. 


O') 


< ||w; K || J2igz + < BjC, 3 , since ||it) rt || < 1. Note that B n < B n+ 1 for n E N from 

oo 

the recurrence relation of the Bell number involving binomial coefficients [36]. This together with the 
inequality above yields the conclusion. □ 


In the next lemma, we study a property of the derivatives of ip K defined as in (14.4|) . 


Lemma 5.3 If f is of Type(//, oo, {0}) for some g E (—1,1) and g E C°°(I ) satisfies Assumption \5A\ 
then there exist two positive constants c and (, > 1 such that for all n > 1, n E N and x E (0,1], 


M / \ 
Vh (x) 


< c2 n B n ( n \rX fJ '~ 


Proof: Applying the Leibniz formula (|4.12|) to yields g^k\x) = YljeZn (/ (X r x))^ vf™ ^\x) 
for x E (0,1] and n E N. By the assumption on /, there exists a positive constant c such that for all 
k > 1, n E N and x E (0,1], (/ (A r x))^ < cXfx^~ n . Using Lemma (15.21) with the inequality above, 
e that there exist two positive constants c and £ > 1 such that for all k > 1, n E N and 

Cl. The desired estimate of this lemma is then obtained from 


we concluc 
x € (0,1], 
this inequality wit 


M / \ 
H>k (x) 


< cB n ( n XrX tJ ‘~ n 

i m . 




□ 


We need two lemmas regarding the parameters rrij for j E Z^ 1 . 

Lemma 5.4 Let s € N with s > 1 and e > 1. If j, l € Z^ 1; j 7 ^ l, then mj 7 ^ m;. Moreover, there holds 
(1 + m)C? - 1) < 2 mj - 2 /or j E 

Proof: We prove the first result by contradiction. Assume to the contrary that there exist j, l E Zj/j with 
j 7 ^ l such that mj = m/. Without loss of generality, we assume j > l. Then, there exists a positive integer 
k such that j = l + k. We then observe that (l + k)e = je < rrij = mi < le + 1, which leads to e < fc _1 < 1. 
This contradicts the assumption e > 1 and thus rrij 7 ^ mi. 

The second inequality follows from the assumptions that — 1 < g < 1 and e > 1, from which we 
conclude that (1 + g)(J — 1) < 2(j — 1 ) < 2 je — 2 < 2 mj — 2 . □ 


Lemma 5.5 If 7 E (0, 1), e > 1 and ( > 1, f/ien f/iere exists a positive constant c such that for all s E N 
with s > 1 and j E Z+_ 1; B 2 mj+i^ 1 +M ^ _ 1 ) + ( 2 m 3 + 1 )/( 2 rrij + 1)! < c, where v := C/ 7 . 

Proof: The proof is similar to that of Lemma 4.1 in [3]. Using the second result of Lemma [5.41 and the fact 
C > 7 , we obtain that z/ 1 +0(i- 1 )+( 2 m j+ 1 ) < Applying inequality ([B.7[i with n := 2rrij + 1, 

we find that 1/(2 mj + 1)! < l/\/2vre (e/(2 mj + l)) 2 m 3 +3 / 2 . Using the bound of Bell number B n < 
(0.792n/ln (n + l)) n with n := 2rrij+l, we have that B 2 mj +i < (0.792(2mj + l)/ln (2m 3 + 2 )) 2m - 7+1 . Com¬ 
bining these three inequalities yields i? 2 mj+i^ 1 + ^^' - 1 ^ 2 mj ’ + 1 V( 2 m .j + 1 )! < d^ inj+1 /(^/2-K(2mj + l)i/ 3 ), 
where dj := 0.792iAe/ In {2mj + 2). It suffices to prove that there exists a positive constant c such 
that for all s E N with s > 1 and j E Z/“_i, d ? J < c. According to the definition of dj and mj, 
we observe that dj < 0.792z/ 2 e/ln(2je + 2). It follows for j E N satisfying j > (2e)~ 1 (e°' 792l ' 2e — 2) 
that < 1. On the other hand, for j E N satisfying j < (2e) _ 1 (e°' 792l/2e — 2), we have that 

d^ rrij +1 < max|d ^ mj+1 : j E N, j < (2e) _ 1 (e°’ 792l/2e — 2)|, which is a constant. This proves the desired 
inequality. □ 

We are now ready to estimate the error £^[f,g] and A f (Q^[/, g]) of the quadrature rule Q^[f,g\. 

Theorem 5.6 Let 7 E (0,1) and e > 1. If f is of Type(/x, 00 , {0}) /or some /i E (—1,1) and g E C°°(I) 
satisfies Assumption 15.11 f/ien f/iere exists a positive constant c such that for all k > 1 and integers s > 1 

£ 7 s [/,<7] < c«; ( ( r 1 ) +M)/(r+1) 7 (1+ ' i)(s - 1) - 1 , 


(5.2) 

















and Af (Q® [/,g]) < |V|(s 2 — s). In particular, if g(x) = x for x € I, then the upper bound in m reduces 
to ck -i-m 7 (i+^)(^-i)-i. 


Proof: We prove (15.21) by estimating £q and Sj[<p K \ := 


X [ Xj ,X j+1 )^ 


r\[ x j i x j+ 1 ] 
W,rrij 


YPk] , for j <E Z+ 1? 


and 


then summing them over j. 

We first consider £o [¥>«]■ By applying Lemma [4.41 with g := A r and b := x\ = 7 s , there exists a 
positive constant c such that for all k > 1 and s € N with s > 1 


£o YPk]< [ \f (X r x)\dx < cXff'j^^ 3 x) . 
Jo 


For j € , by the error bound of the Gauss-Legendre quadrature, there exists a constant € [xj,x J+ \] 


such that £j[<p K } < /i • mj /( 2 2m ? (2mj + 1 )!) 


A 2m,) te: 


where hj := Xj+i — 27 . Applying Lemma [5.31 with 


> Xj for j € ZJ 1 1 , there exist two positive constants c and £ > 1 such that for all k > 1 , s € N with 


s > 1 and j € Z^_ 1 , 


A 2m,) &: 


< c2 2m i 


Combining these two inequalities yields two 


positive constants c and ( > 1 such that for all k > 1, s € N with s > 1 and j € Zj "_ 1 


Ro r 2mu2™,+l/-2m, R 

£ , w < cAf MA_2>— < c B2m ' +1 " 


( 2 mj + 1 )! 


{2m,j + 1 )! 


-_\Af 7 ( 1 +M)(«- 1 )(l _ 7 )2mj+l 


From Lemma 15.51 we conclude that there exists a positive constant c such that for all k > 1, s € N with 
s > 1 and j € Z^_ 1 , < cAr 7 ^ 1+At ^ s_1 ^(l — 7) 2m,+1 . Summing up the bound of errors £j[<p K \ over 

j € Z s _i, we conclude that there exists a positive constant c such that for all k > 1 and s € N with s > 1, 
£ 8 [/, g\ < cA £ + V 1+M)(s_1) (l + S 7 ez+ (■*■ — 7 ) 2 m j+ 1 j. According to my 7 ^ mi for j,l € Zjl x with j / l 

proved in Lemma 15.41 we observe that {2 my + 1 : j € Zj~_[} U {0} C Z := {0,1,2,...}. It follows that 
1 + SjGZ+ (1 — 7 ) 2 m J +1 < Yljs z(! — 7) J ” = 7 _1 - Substituting this result and the definition (14.21) of A r into 
the inequality above yields the estimate (15.21) . 

The bound of the number of functional evaluations used in quadrature (15.11) may be obtained directly. 
When g(x) = x, we substitute r = 0 and k ct (o) = k into estimate (15.21) to yield the special result. □ 


We next develop a quadrature rule for the oscillatory integrals (14.51) . This is done by choosing variable 
numbers of quadrature nodes in the subintervals of A. Specifically, for n € N and for the partition n K(7 ( r ) 
of A with nodes defined by (14.161) . we let rrij := n + \(n + 1 — j)( 1 — a)] and Nj := \qf\ for j € Z+. For 
each j € Z+, we use [/) ( j\ t0 approximate the integral x\f' J ^ uX " J \f,g\. Integral Xjf[f,g] defined by 

(14.51) is then approximated by the quadrature formula 


Q*,n[/)S] : = 


Ec 

je z+ 


[Xj-l,Xj] 




(5.3) 


The next lemma provide an estimate of the errors £j [/, g] 



[/> 9} - Q 


[Xj~l,Xj]r j. n 

Nj,K,mj J 


for j € Z+. 


Lemma 5.7 If 0 is of Type(a, 00 , {0}) for some a € (—1,1) and g satisfies Assumption HJ], then there 
exists a positive constant c such that for all k > 1, n € N and j € Z+ 


£?[/, 5 ] < c 


iff(i)i~ Q mr - 1 

( 1 rrij - l)\K 2 N™ j l 


(<*-l)(r+l) 

®J -1 


(A*(r)) m ' • 


This lemma may be proved in the same way as Lemma 14.61 with m = mj. We next estimate the error 
£ K,n [/> ll] := |Z« [/,g] - To this end, we establish two technical lemmas. 


Lemma 5.8 For all n £ N with n > 1, z + 1 /(m 3 - — 1) < 1. 















Proof: The proof is similar to that of Lemma 13.31 For n € N with n > 1, we need to estimate the lower 

bound of the set {rrij — 1 : j € Z+}. By a < 1, we have that rrij — 1 > n for j € Z+. It follows that 

l/( m i — 1) ^ V n = 1- ^ 

For simple notation, we let r n (r) := ~ l). 

Lemma 5.9 There exists a positive constant c such that for all At > 1, j G Z+ and n € N satisfying 

r n (r) < (n- l)/e, (5.4) 

f/iere holds (T n (r)) mj ~ n /(rrij — 2)! < c(n — l) -1 / 2 . Moreover, a t*^ < 2(n— l) 1 / 2 . 

Proof: By a < 1 we see that nij — 2 > n — 1 for j € Z+. Condition (15.41) implies n > 2 since Ti(r) > 0 for 
k > 1. It follows for n € N satisfying (15.41) with r n (r) < 1 that (r n (r)) mj_n /(rrij — 2)! < 1 /{rrij — 2)! < 
(n— 1) _1//2 . On the other hand, for n £ N satisfying (15.41) with r n (r) > 1, we have ( T n (r)) mj ~ n /(mj — 2)! < 
( T n{r)) m:> ~ 2 /(m,j — 2)!. According to inequality (13.71) with n := mj — 2, there exists a positive constant 
c such that for all n € N with n > 1 and j € Z+, l/{mj — 2)! < c(mj — 2) -1 / 2 ( e/(m,j — 2)) mj ~ 2 < 

c(n — l ) -1 / 2 (e/(n — l)) mj ~ 2 . This together with the inequality above ensures that there exists a positive 

constant c such that for all At > 1, j € Z+ and n € N satisfying (15.41) with r n (r) > 1, ( T n (r)) mj ~ n / (rrij— 2)! < 
c(n — l ) -1 / 2 (er n (r)/(n — l ))” 1 -? -2 < c(n — l)” 1 / 2 . This concludes the desired result. 

We next prove the second inequality. From condition (15.41) we observe that — — T n( r ) < 

(n — l)/e. Using the inequality ((n — l)/e ) 1 ^ 2 + 1 < 2(n — l ) 1 / 2 for n > 1 yields that for n G N satisfying 
(15.41) with > 2, At^” < ((n— l)/e ) 1//2 + 1 < 2(n — l) 1 / 2 . Moreover, it follows for n € N satisfying 

(|5.4D with < 2, At^jf < 2(n — l) 1 / 2 . Combining these two results we obtain the second conclusion. 

□ 

We are now ready to estimate the error £^ n [f,g] and A f (Q« n[/, <?]) • Let p n (r) := 1 — A X J n , where A r 
is defined by (14.21) . We first consider the case r = 0. For simple notation, let := r n (0) and := p n ( 0). 


Theorem 5.10 7/0 is of Type(a, oo, {0}) for some a € (—1,1) and a satisfies Assumvtion \-L.l\ with r = 0. 
then there exists a positive constant c such that for all k > 1 and n € N satisfying m, 


£ K,n\fi9\<c\g( 1)1 Q o- 0 <5o 2 (n — 1) 1/2 a e 2 At ff0 (p°)". 

For n € N, i/iene holds the estimate N (Q^ n [/, < 7 ]) < [oo/^ol (2n 2 + |"1 — a\ (n 2 + n))/2 + 1. 


(5.5) 


Proof: We prove (15.51) by estimating £j [/, g] and then summing them over j G Z+. By employing Lemmas 
I4.7l and l5.7l with r = 0, and the choice of Nj , there exists a positive constant c such that for all At > 1, j G Z+ 

and n G N, [/, a/] < - iJL—Applying = A 0 ] ^ n p[{ with the definition of Xj- 1 , 

(mj — 1)!ac 

we observe that x^Z\ (/3n) mj < ” Aq 1 (p°) n . Combining these two inequalities yields a positive 

\g(l)\~ a aoS%~ 2 (t°)" 


constant c such that for all At > 1, j G Z+ and n G N, £j[/, g] < c- 


(Ao 1 K) n 


(mj — l)At 2 (mj — 2)! 

This together with the first result of Lemma 15.91 with r = 0 ensures that there exists a positive constant 
c such that for all At > 1, j G Z+ and n G N that satisfy (15.41) . £j[f,g] < c(mj — 1 ) _1 |< 7 (l)p tt ctqSq~ 2 (n — 
l)^ 1 / 2 At _2 AQ 1 (pn) n ■ Summing up the inequality above over j we obtain that there exists a positive 

constant c such that for all At > 1 and n € N satisfying m , 

s £,nlf>9\ < c\g(l)\~ a a 0 5f- 2 (n - l)~ 1/2 h~ 2 \q 1 (p° n ) n ^ (mj - l) -1 . 


This together with Lemma 15.81 and A 0 1 = 


K ao leads to the estimate (15.51) . 




According to formula (15.31) and Lemma 14771 with r = 0, we obtain that for n G N, Af (Q^ n [/, g]) < 
{ r°oA>l m j + i} - (n- i) < f<To/<5 0 ] (2n 2 + [1 - a](n 2 + n ))/2 + 1. □ 

Following Theorem 15.101 we have a special result for the case g{x) = x for x G I. 

Corollary 5.11 If f is of Type(/x, oo, {0}) for some /x € (—1,1), f/ien f/iere exists a positive constant c 
such that for all k > 1 and n G N that satisfy (15. 4|) with r = 0, £*„[/, 3 ] < c(n - 1) 1 / 2 k 1 e((, where 
e n ■= 1 — kT 1 /”. For n G N, t/iere holds the estimate Af g ]) < (2n 2 + |~1 — /x] (n 2 + n))/2 + 1. 

Next, we consider the case r > 0. 


Theorem 5.12 7/0 is o/Type(a, oo, {0}) for some a G (—1,1) and g satisfies Assumption[ XT] with r > 0, 
then there exists a positive constant c such that for all k > 1 and n G N satisfying m, 

£^ n [f,g] <c(r + l)\g{l)\~ a a{r){S{r)) a - 2 K- 2 Kl / ( ^ +1 \p n (r)) n . (5.6) 


For n G N, t/iene holds Af (Q^ n [/, <y]) < 


(r + !) 


r/(n(r+l)) 

'MO 


cr(r)/(5(r) (2n 2 + |~1 


— a] (n 2 + n))/2 + 1. 


Proof: The proof is handled in the same way as that of Theorem 15.101 Using /3 n {r) = A r 1 ^ n pn{f) and 
the definition of Xj- 1 , we obtain that {/3 n {r)) m i < {T n {r)) mj ~ n Xf 1 {p n {r)) n , for j G Z+. By 

employing Lemmas 14.71 and 15.71 the choice of Nj and the inequality above, there exists a positive constant 
c such that for all k > 1, j G Z+ and n G N, 


(r + l)|g(l)| “ a{r)(6{r)) a 2 K 1 ( n . ( T (r)) mj ^ n 

+/.«] < 4 — „ , : — ” „„ Klpnir))’ 


{mj — l)n 2 


{mj - 2 )! 


Using the first result of Lemma 15.91 we obtain that there exists a positive constant c such that for all k > 1, 
j G Z+ and n G N satisfying (|5.4|> . 

£?[/,#] < c(r + 1) |5(l)r Q a(r)( ( 5(r))"” 2 K^ ) (n - l)” 1/2 K" 2 A“ 1 (p„(r)) n /(m i - 1). 

This together with the second result of Lemma 15.91 and A^ 1 = ensures that there exists a positive 

constant c such that for all k > 1, j G Z+ and n G N satisfying dSZD, 

< c(r + 1 )| 5 (l)r a a(r)(5(r)) Q - 2 «;- 2 ^ +1) (p n (r)) n /(mj - 1). 

Summing up the inequality above over j G Z+ and applying Lemma 15.81 we obtain the estimate (15.61) . 

It remains to estimate Af {Q£ n [/, 5 ]) • According to formula ()5 .31) and Lemma [4.7l we have for n G N that 

^{Qin[f,g]) < E ieZ + { fal mj + l}-(n-l) < (r + l)K^ (r+1)) <T(r)/5(r) (2n 2 +[l-a](n 2 +n))/2+l. 
This yields the last conclusion. □ 


Note that the decay of the bound of £^ n [f,g\ given in Theorems 15.101 and 15.121 is faster than the 

exponential decay with the base p no { r )> where no := min{n : n G N,n satisfying condition (|5.4p l. For 

k > 1, we conclude that 0 < /9 nn (r) < 1 and lim {pn{r)) n /{p nn { r )) n = 0. 

n—>• 00 

The oscillatory integral (12.11) may be computed by Q'lfn[f- g] := Q®[/, 5 ] + Q« n [/, S']- We next es¬ 
timate the error £ff [/, g\ of the quadrature formula Qnji[f,g] under reasonable hypotheses. Note that 
—n In (l — 7 1+/i ) > In n CT ( r j for n G N implies that p n {r) < 7 1+/ L If r > 0 and n = s, then there exists a 
positive constant c such that for all k > 1 and n G N that satisfy (15.41) and —nln (1 — 7 1+At ) > In /+(,.), 
£f’f [/, g\ < If r = 0 and n = s, we consider two cases. In the case —1 < p < 0, 

there exists a positive constant c such that for all k > 1 and n G N satisfying m and —nln (l — 7 1+/i ) > 
In , £f{n [/, g] < . In the case 0 < /x < 1, there exists a positive constant c such that 

for all k > 1 and n G N satisfying (15.4(1 and —nln (l — 7 1+At ) > ln/c CT0 , £2’,n[f,g] < CK~ 1 'y( 1+ ^( n ~ 1 ' > . 















6 Numerical Experiments 


In this section, we present numerical results of seven computational experiments to verify the approxi¬ 
mation accuracy and computational efficiency of the proposed composite moment-free Filon-type (CMF) 
quadrature formulas, among which the quadrature formulas with a polynomial order of convergence will be 
denoted by CMFP and the quadrature formulas with an exponential order of convergence will be denoted 
by CMFE. The numerical experiments to be presented are divided into two groups. The first three exam¬ 
ples are about the CMFP formulas and the last four examples concern with the CMFE formulas. We also 
compare the computational performance of the proposed quadrature formulas with that of the quadrature 
rules proposed in mm- The numerical results presented below are all obtained by using Matlab in a 
modest desktop (a Core 2 Quad with 4Gb of Ram memory). 

In the first example, we test the CMFP formula (I3.4|) for computing the oscillator integral with smooth 
integrand and linear oscillator. 


Example 6.1 The purpose of this example is to confirm the theoretical estimate presented in Proposition 
o for the CMFP formula (|3.4I) which approximates the integral Z K [/, g\. In this example, we consider 
the functions f(x) := e x and g(x ) := x for x € I. The exact value of the corresponding integral can be 
computed exactly, and thus, the errors presented below are computed using the exact value. 

Table 6.1: Numerical results of CMFP for f(x) := e x and g(x) := x 


K, 

n = 5 

n = 15 

n = 20 

n = 30 

RE 

AT 

RE 

N 

RE 

JV 

RE 

AT 

l(F 

4.68e-7 

21 

8.93e-8 

61 

4.13e-8 

81 

1.15e-8 

121 

10 3 

2.43e-7 

21 

2.21e-9 

61 

1.24e-8 

81 

1.95e-9 

121 

10 4 

4.62e-8 

21 

2.99e-9 

61 

1.77e-9 

81 

1.62e-10 

121 

10 s 

5.72e-9 

21 

6.38e-10 

61 

3.47e-10 

81 

1.65e-ll 

121 

10 6 

1.36e-10 

21 

2.62e-10 

61 

1.26e-10 

81 

2.67e-12 

121 

10 7 

4.54e-ll 

21 

1.51e-ll 

61 

5.93e-12 

81 

2.79e-12 

121 



A: RE of Q K , n Af,9] 

Figure 6.1: RE for f(x) : = 


B: RE of Q K , nA [f,g] 
e x and g(x) := x 


We present in Table lfTTI the relative error (RE) of Q K n ,A\fi s\ an d th e number AA of functional evaluations 
used in the formula for different values of k and different choices of n. We plot in Figure l6Tl the RE values 
of the quadrature formula for (A) n = 10 3 and for (B) k = 10 4 . Figure l6Tl and Table I7T7T1 confirm that for 
a fixed k the approximation accuracy of the quadrature formula increases as n grows and for a fixed n it 
increases as k grows. 

We verify in Figure lR2l the order in 1/k of the error £ K , 8 , 4 [/,<?]• We plot the error £ K ,8,4 [fid] scaled by 
k 2 in (A) and that scaled by k 3 in (B), for n = 8 with k changing from 10 2 to 10 3 . Comparing (A) and 
(B) of Figure 1(01 we observe that the decay of the error £ K , 8 , 4 [/, 5 ] is faster than 0(k~ 2 ), but is slower 
than 0(k~ 4 ). 































A: K 2 S Ki8A [f,g] 


B: k 3 £, 


«,8,4 


[/, 5] 


Figure 6.2: Error £ Kj g, 4 [/,<?] scaled by k 2 and k 6 for f(x) := e x and g(x) := x 


We next compare the performance of the proposed quadrature formulas with that of the existing 
formulas described in El d- To this end, we recall the quadrature formulas of E1E|. Following El E|, 
by f3 £ (—1,1) we denote the index of singularity of /, N € N denote the degree of the polynomial 
interpolant, q > 0 denote the grading parameter and M £ N denote the number of the subintervals. The 
Filon-Clenshaw-Curtis (FCC) rule proposed in |6| for the integral xl g] with a smooth integrand 

/ and a linear oscillator g was designed by approximating / by its Lagrange interpolation Qn(I) of 
degree N at the Clenshaw-Curtis points tj := cos (jn/N) for j £ TLjq. The FCC rule proposed in [6] 
for the integral m with a smooth integrand and a linear oscillator (that is, g(x) = x) is given by 
1 !°jv (/) := e 1K / 2 zj.y 2 ’^[QAf(/),5']/2, where f(x) := f((x + l)/2) for x £ [—1,1]. When / has an integrable 
singularity at the origin, a composite Filon-Clenshaw-Curtis (CFCC) quadrature rule was proposed in [TJ. 
Specifically, a positive integer M is chosen, q that satisfies the condition q > (7V + 1 — r)/(l + /3 — r) for some 
0 < r < 1+ /3 is picked, and a partition of the interval I is formed with nodes Xj := M~ q j q for j £ Zm- On 
all (but the first) subintervals of the partition, the FCC rule is applied. The treatment of the first subinterval 
depends on the value of the parameter /3. For (3 £ (—1,0], the integral x)// 0,Xl ^ [/, g\ is approximated by zero. 
For /3 £ (0,1) with x \k > 1, the linear function interpolating the two pints (xo,/(xo)) and (xi,f(x\)) is 
used to approximate the function / on [xo,a:i]. While for (3 £ (0,1) with x\K < 1, the integrand of (12.11) on 
[xq,xi\ is approximated by the linear function interpolating (xo,/(xo)e 1K5 Oo)) and (aq,/(xi)e 1K5 Oi)). The 
FCC rule with same number N + 1 of points is used to calculate the integrals defined on the subintervals 

[xj,Xj. |_i] for j £ Let ij °’ ^(/) denote the approximation of the integral , g\. The CFCC 

rule is formed by X [ ^ ] Mq (f) := Tj^ ll (/) + £ jeZ +_ i g\, where hj := (x j+1 -Xj)/2, 

dj := (xj + 1 +Xj)/ 2 and fj(x) := f(hjX + dj) with x € [—1,1] for j £ When g is a strictly increasing 

nonlinear function and has a single stationary point at zero, we calculate the integral (12.31) instead by a 
change of variables g(x) —>• x and mapping [g(0),g(l)] onto I. The quadrature formulas proposed in El Ei 
are implemented by using the Matlab code available in [8j- 

We now discuss the computation of the errors to be presented in Examples 16.2116.3116.51 and 16.61 The 
true values of the integrals considered in these examples cannot be computed exactly. Therefore, we need to 
find their appropriate substitutions. This is done by the means of the Gamma function. To this end, we let 
Z K [0] := /q 1 lnxe 1K:r dx and X K [a\ := x a e 1KX dx, for a € (—1, 0) U (0,1). According to [IT] (see also |7|), we 
have that X K [0] = — (7 t/2 + L + ilnK) /k+0{1/k 2 ) andX K [a] = (T(l + a)e 1,r ( 1+ “)(ifi:) - ( 1+Q ) + e 1K /(in)) (1+ 
0(1/k)), for a £ (—1,0) U (0,1), where T denotes the Gamma function, arg(i«:) = 7r/2, and ? denotes the 
Euler-Macheroti constant having an approximate value 0.5772. We shall use 1/ [0] := — (vr/2 + i? + i In k) /k 
and X“[a] := T(1 + a)e 1,r ( 1+ ")(iK) - ( 1+Q ) +e 1K /(i«:), for a € (—1,0) U (0,1), respectively, to approximate the 
exact values of X K [0] and X K [a \, and substitute them when we compute errors of quadrature rules. 
































Table 6.2: Comparison of CMFP and CFCC for fi(x) := x 1//2 and g(x) := x 


K, 

CMFP 

CFCC 

n = 5 

n= 10 

M = 

7 

M = 14 

RE 


RE 

AT 

RE 

JV 

RE 

JV 

l(F 

5.03e-3 

37 

5.10e-3 

77 

5.09e-3 

38 

5.09e-3 

80 

10 3 

5.80e-4 

37 

5.12e-4 

77 

5.65e-4 

38 

5.14e-4 

80 

10 4 

8.40e-5 

37 

5.41e-5 

77 

1.52e-4 

38 

5.10e-5 

80 

10 s 

8.58e-5 

37 

6.94e-6 

77 

7.65e-5 

38 

6.33e-6 

80 

10 6 

3.26e-5 

37 

5.24e-6 

77 

2.17e-4 

38 

3.06e-6 

80 

10 7 

3.02e-5 

37 

2.55e-6 

77 

3.84e-4 

38 

8.58e-6 

80 


In the next two examples, we test the efficiency of the quadrature rule proposed in Section 4 for 
calculating the oscillatory integrals with singular / and with/without a stationary point of g. 

Example 6.2 This example is designed to verify the theoretical estimates given in Theorems 14.51 and 14.81 
for the CMFP formula [/, g\ , which approximates the integral Z K [f,g\. We consider three functions 

/i(x) := x 1 / 2 , f 2 (x) ■= hi x. / 3 (x) := x ” 1 / 2 and g(x) := x for x € I. When we compute the errors of the 
quadrature formulas, the true values of the integrals I K [fj,g\ are computed by using for j € Z$, 

respectively, where a.\ := 1 / 2 , a% := 0 and 0:3 := — 1 / 2 . 


Table 6.3: Comparison of CMFP and CFCC for / 2 (x) := lnx and g(x) := x 


AC 

CMFP 

CFCC 

n = 5 

n = 10 

M = 

7 

M = 14 

RE 

J\ 

RE 

IT 

RE 

Af 

RE 

J\ 


1.75e-3 

37 

1.86e-3 

77 

4.15e-3 

37 

1.86e-3 

79 

10 3 

1.72e-3 

37 

1.19e-4 

77 

8.87e-3 

37 

1.93e-4 

79 

10 4 

2.90e-3 

37 

1.15e-4 

77 

1.59e-2 

37 

1.65e-4 

79 

10 5 

6.23e-3 

37 

1.42e-4 

77 

5.14e-3 

37 

2.45e-4 

79 

10 6 

1.08e-2 

37 

9.16e-4 

77 

4.76e-2 

37 

3.99e-4 

79 

10 7 

1.60e-2 

37 

1.33e-3 

77 

1.53e-2 

37 

1.56e-3 

79 


Numerical results of this example are presented in Tables lfT2ll6.3l and l6. 41 We compare the RE values of 
the formula 4 [/ 1 , g\ and the CFCC formula 8 (/i) in Table HOI those of the formula Q™’n d\ 

and the CFCC formula 12 (h) i n Table [6T3l and those of the formula 4 [/ 3 , 5 ] and the CFCC 

formula xM, i 6 (/ 3 ) in Table [6741 These numerical results show that overall the CMFP formula has higher 
approximation accuracy than the CFCC rule. 


Table 6.4: Comparison of CMFP and CFCC for / 3 (x) := x 1//2 and g(x) := x 


AC 

CMFP 

CFCC 

n = 5 

n = 10 

M = 

7 

M = 14 

RE 


RE 

AT 

RE 

A f 

RE 

AT 

10^ 

2.89e-2 

37 

1.25e-3 

77 

1.18e-2 

37 

4.07e-4 

79 

10 3 

2.50e-2 

37 

9.23e-4 

77 

1 .Ole-2 

37 

6.62e-4 

79 

10 4 

3.14e-2 

37 

7.49e-4 

77 

2.57e-2 

37 

1.05e-3 

79 

10 s 

3.38e-2 

37 

5.54e-4 

77 

4.49e-l 

37 

4.82e-3 

79 

10 6 

8.55e-2 

37 

2.22e-3 

77 

1 .12e-l 

37 

1.70e-2 

79 

10 7 

1 .12e-l 

37 

5.93e-3 

77 

2.06 

37 

4.02e-2 

79 


Example 6.3 This example is to verify the estimates proved in Theorems 14.51 and 14.101 for the CMFP 
formula Q^m [/, g], which approximates X K [f, g\. In this example, we consider the functions f(x) := x 172 
and g(x) := x 2 for x € I. When we compute the errors of the quadrature formulas, the true value of this 
integral is computed by using /“[—3/4]/2. 














































A: RE of Q”’f 4 [/,<?] 


B: RE of Q"’ 4 [/,<?] 


Figure 6.3: RE for f(x ) := x 1 / 2 and g(x) := x 2 




A: AA(Q^ !4 [/, 5 ]) B: RE of Q n K A nA [f,g\ 

Figure 6.4: Af and RE for /(x) := x -1 / 2 and g(x) := x 2 


Numerical results of this example are presented in Figure [6T3V16.61 and Table [6751 We plot in Figure 16431 
the RE values of the formula Q”’ 4 4 [/, g] for (A) k = 10 3 and for (B) k = 10 4 . We present in Figure EH 
the number A f of functional evaluations used in the formula 2^5 4 [/, d\ (A) and the RE values (B), and 
those of the formula 2^30,4 [/i #] i n Figure ESI We plot in Figure E6l the RE values obtained from the 
CMFP and CFCC formulas by using the same number of functional evaluations for (A) k = 10, for (B) 
k = 10 5 and for (C) k = 10 6 . The RE values of the formula Q”d_ 4 [/, g\ and those of the CFCC formula 

-Z>,° 4 M 21 ((/ /s') °ff^) are reported in Table ESI 



A: 


B: RE of Q n K f nA [f,9] 


Figure 6.5: A( and RE for /(x) := x 4 / 2 and g(x) := x 2 


From these numerical results, we have the following conclusions. For a fixed k the approximation 
accuracy of the CMFP formula increases as n grows. The CMFP formula achieves an accuracy of an order 
higher than the CFCC formula for large k and is effective and convenient for implementation. 




















































Table 6.5: Comparison of CMFP and CFCC for f(x) := x 4 / 2 and g(x) := x 2 


K, 

CMFP 

CFCC 

n = 30 

n = 50 

M = 200 

M = 

400 

RE 

AJ 

RE 

AT 

RE 

AT 

RE 

Al 

l(F 

6.62e-5 

477 

6.59e-5 

797 

6.60e-5 

797 

6.61e-5 

1597 

10 3 

1.47e-6 

477 

1.17e-6 

797 

8.81e-7 

797 

9.50e-7 

1597 

10 4 

4.33e-7 

597 

1 .21e-8 

797 

6.016-7 

797 

3.42e-7 

1597 

10 s 

4.52e-7 

597 

8.95e-9 

797 

4.57e-7 

797 

3.92e-7 

1597 

10 6 

4.55e-7 

597 

9.31e-9 

997 

2.91e-7 

797 

3.07e-7 

1597 

10 7 

4.45e-7 

597 

8.88e-9 

997 

5.64e-7 

797 

3.72e-7 

1597 





number of evaluations number of evaluations number of evaluations 


A: K = 10 4 B: k = 10 5 C: k = 10 6 

Figure 6.6: Comparison of CMFP and CFCC for f(x) := £ -1 / 2 and g(x) := x 2 


The numerical results of these three examples indicate that the CMFP formulas are stable, easy to im¬ 
plement and computationally inexpensive. For irregular oscillator, the performance of the CMFP formulas 
is superior to that of the FCC and CFCC formulas. The advantage of the proposed formulas is especially 
obvious for the oscillator with stationary points. We next show the numerical results calculated by the 
CMFE formulas. 

In Example 16.41 we test the CMFE formula (13.51) for computing the oscillator integral considered in 
Example 16.11 and compare with the CMFP formula (13.41) . 

Example 6.4 This example is to confirm the estimate shown in Theorem 13.51 for the CMFE formula 
(13. 5 p which approximates the integral I K [f,g]. We consider in this example the functions f(x) := e x and 
g(x) := x for x €E I, the same as those in Example 16.11 

Table 6.6: Numerical results of CMFE for f(x) := e x and g(x) := x 


K 

n = 4 

n = 5 

RE 

AT 

RE 

AT 

10 3 

9.65e-14 

26 

2.24e-13 

47 

10 3 

2.17e-15 

26 

6.39e-15 

47 

10 4 

7.36e-17 

26 

1.62e-15 

47 

10 s 

1.14&-17 

26 

5.70e-17 

47 

10 6 

1.30e-16 

26 

2.61e-16 

47 

10 7 

1.50e-16 

26 

2.32e-16 

47 


Numerical results of this example are reported in Table I7TTTT1 and Figure [6771 We list in Table ltT6l the RE 
values of the formula Q KjU [f, <?]• We depict the RE values of the formula Q Kt n[f,g] in Figure [R7l (A) and 
the error £ KtU [f,g\ scaled by K n+1 in Figure RTTl (B), for n = 4 with k changing from 100 to 5000. From 
Table [TTTTTl and Figure l6?7l we observe that the approximation accuracy of the formula Q re ,n[/> s\ increases as 
k grows for a fixed n and the asymptotic order of convergence is 0{n~ n ^ 1 ) for n = 4. It concurs with the 
theoretical estimate. Comparing Figure I67F1 (B) with Figure [6721 in Example 16.11 we can obtain an order of 
convergence faster than 0 (k~ 3 ) by using variable number of quadrature nodes in the subintervals. 




















































A: RE of Q K , n [f,g] B: K n+1 £^ n [f, g] 

Figure 6.7: RE and error S K ,n[f,g] scaled by K n+1 of Example 16.41 


In the following two examples, we validate the efficiency of the quadrature rule proposed in Section 5 
for calculating the oscillatory integrals with a stationary point, without/with singular /. 

Example 6.5 This example is to verify the estimates established in Theorems 15.61 and 15.121 for the CMFE 
formula QkJi [f, g] , which approximates the integral I K [/, g]. We consider in this example the functions 
f(x) : = 1 and g{x) := x 3 for x € I. When we compute the errors of the quadrature formulas, the true 
value of this integral is computed by using /“[—2/3]/3. 


Table 6.7: Comparison of CMFE and CFCC for f(x) := 1 and g(x) := x 3 


K 

CMFE 

CFCC 

n = 

3 

n = 

4 

M = 

400 

M = 

800 

RE 

TV 

RE 

M 

RE 

N 

RE 

Al 

ted 

1.17e-4 

307 

1.17e-4 

365 

1.18e-4 

1597 

1.17e-4 

3197 

to 3 

2.49e-6 

407 

2.48e-6 

467 

2.62e-6 

1597 

2.43e-6 

3197 

to 4 

5.36e-8 

607 

5.36e-8 

603 

5.07e-7 

1597 

3.04e-7 

3197 

to 5 

1.08e-9 

907 

9.38e-10 

841 

5.82e-7 

1597 

3.76e-7 

3197 

10 6 

9.45e-ll 

1427 

4.00e-10 

1156 

7.22e-7 

1597 

4.04e-7 

3197 

to 7 

7.13e-ll 

2287 

7.22e-10 

1657 

6.10e-7 

1597 

4.09e-7 

3197 


Numerical results of this example are presented in Table 16771 and Figure [6781 We list in Table 16771 the 
RE values of the formula QLfn 02 [f, g\ and those of the CFCC formula 16 (( f/g ') o In (A) of 

Figure [6781 we depict the RE values of the formula Q'^ 3 ' 02 [/, g\ , which are listed in Row 2 of Table [6771 and 

the RE values obtained from the CFCC formula zj . 0 ^ 1 ](j {[f/g') o^ 1 )), where the numbers AT = 309, 

409, 609, 909, 1429, 2289 of functional evaluations used in the formula are comparable to those used in 
the CMFE formula recorded in Row 3 of Table [673 In (B) of Figure [6781 we compare Af{Q 7 f° 3 02 [f, g]) with 
(r + 1 )n 2 «: 1//n with r = 2 and n = 3. 



A: RE of Q 7 S y2 [f\g] 

Figure 6 . 8 : RE and A f for f{x) 



B: Af(Q 7 ’° 3 02 [f,g]) 
:= 1 and g{x) := x 3 

















































Table 6.8: Comparison of CMFE and CFCC for f{x) := x 3 / 2 and g(x) := x 2 



CMFE 

CFCC 

n = 

3 

n = 

4 

M = 100 

M = 

300 

RE 

M 

RE 

Af 

RE 

JV 

RE 

AJ 

10 * 

6.59e-5 

502 

6.59e-5 

537 

6.48e-5 

397 

6.58e-5 

1197 

10 3 

1.19e-6 

544 

1.16e-6 

572 

3.29e-6 

397 

1.47e-6 

1197 

10 4 

6.23e-8 

607 

1.75e-8 

642 

6.56e-6 

397 

4.52e-7 

1197 

10 5 

4.50e-8 

691 

2.74e-9 

712 

1.04e-5 

397 

3.28e-7 

1197 

10 6 

1.90e-8 

829 

2.13e-9 

817 

1.77e-5 

397 

4.93e-7 

1197 

10 7 

1.09e-8 

1027 

4.12e-9 

922 

3.15e-5 

397 

4.45e-7 

1197 


We conclude from the results in Table [6771 and Figure [6T8l (A) that for a fixed n the approximation accu¬ 
racy of the formula Q;-'.^ 02 [/, g\ increases as k grows and AT(Q 7 K ’° n 02 [f, g]) increases slow as k grows. Further¬ 
more, from the results presented in Table 1(771 we observe that for n = 3, 4 and large k. quadrature formula 
A f (Q 7 k%° 2 [f, g\) has better approximation accuracy than the CFCC formula 16 ((/ /g') o , which 

uses more number of functional evaluations. In fact, from Figure [6781 (B). we see that the number of func¬ 
tional evaluations used in formula Q^ , .!' 02 [/, g] is significantly less than 27 k 1 / 3 for large k. 

Example 6.6 The purpose of this example is to test the estimates proved in Theorems 15.61 and 15.121 for 
the CMFE formula Q'J,n[/, g]■ which approximates the integral X K [f,g\. In this example, we consider the 
functions fix) := x~ l f 2 and g(x) := x 2 for x € I, the same as those in Example 16.31 



We list in Table 16.81 the RE values of the formula Ql^n ' 02 [f, g\ and those of the CFCC formula 
^I°4M2i (( f/g') ° g I 11 (A) °f Figure [FTTH we depict the RE values of the formula Q^ 2 3°' 02 [/, <?], which 
are listed in Row 2 of Table [6781 and the RE values of the CFCC formula 21 (( f/g') o gU 1 )) in (A) 

of Figure ITITTJl where the numbers Af = 505, 545, 609, 693, 829, 1029 for k = 10 J+1 , j € of functional 
evaluations used in the formula are comparable to those used in the CMFE formula recorded in Row 3 of 
Table lR8l In Figure I(u9l (B), we compare Af{Q 12 f' 02 [f, g]) with (r + 1 )n 2 K 1 f n with r = 1 and n = 3. 

From Figure [6791 (A) and Table 1(781 we conclude that the accuracy of g\ is higher than that 

of the CFCC formula. 

From Examples 16.51 and 16.61 we observe that the CMFE formula proposed in Section 5 has higher order 
of approximation accuracy than the CFCC formula when the integrand has a singularity with the index 
g, < 0 and a stationary point of the order r > 0. Moreover, the CMFE formulas do not have to compute 
g- 1 for the nonlinear oscillator g (which is required by the FCC and CFCC formulas) and as a result, they 
can save an enormous amount of computing time in comparing with the FCC and CFCC formulas. 

In the next example, we shall compare the CPU time spent using the CMFE formulas, the FCC and 
CFCC formulas when computing oscillator integrals with a nonlinear oscillator. The CPU time is monitored 
by using tic and toe of Matlab. The computation of the change of variables g(x) —>• x is carried out using 
fzero of Matlab. 




































Example 6.7 This example compare the CPU time spent when computing the integral I K [f,g] by using 
the CMFE formulas, the FCC and CFCC formulas. We consider the functions fi(x) := 1, f 2 (x) := lnx 
and g{x) := (sin (nx/2) + 2x) /3 for i£l, 


Table 6.9: Comparison of CPU time for fi(x) := 1 and g{x) := (sin (irx/2) + 2x) /3 



CMFE 

FCC 


Time (sec.) 

Approximation 


Time (sec.) 

Approximation 

J\ 

10* 

3.34e-2 

-7.35e-3-(4.66e-3)i 

12 

2.84e-l 

-7.35e-3-(4.66e-3)i 

9 

10 3 

3.35e-2 

1.24e-3-(1.12e-6)i 

12 

2.97e-l 

1.24e-3-(1.12e-6)i 

9 

10 4 

3.31e-2 

-4.59e-5+(2.27e-4)i 

12 

2.85e-l 

-4.59e-5+(2.27e-4)i 

9 

10 5 

3.33e-2 

5.36e-7+(2.34e-5)i 

12 

2.83e-l 

5.36e-7+(2.34e-5)i 

9 

10 6 

3.31e-2 

-5.25e-7-(5.65e-7)i 

12 

2.84e-l 

-5.25e-7-(5.65e-7)i 

9 

10 7 

3.35e-2 

6.31e-8+(2.20e-7)i 

12 

2.85e-l 

6.31e-8+(2.20e-7)i 

9 


Table 6.10: Comparison of CPU time for f 2 (x) := In x and g{x) := (sin (nx/2) + 2x) /3 


K 

CMFE 

CFCC 

Time (sec.) 

Approximation 

AT 

Time (sec.) 

Approximation 

JV 

10* 

4.09e-2 

-1.30e-2-(4.51e-2)i 

88 

3.90e-l 

-1.30e-2-(4.51e-2)i 

73 

10 3 

4.10e-2 

-1.31e-3-(6.43e-3)i 

88 

3.95e-l 

-1.32e-3-(6.43e-3)i 

73 

10 4 

4.10e-2 

-1.32e-4-(8.37e-4)i 

88 

3.96e-l 

-1.31e-4-(8.37e-4)i 

73 

10 5 

4.10e-2 

-1.32e-5-(1.03e-4)i 

88 

3.63e-l 

-1.33e-5-(1.03e-4)i 

73 

10 6 

4.15e-2 

-1.32e-6-(1.22e-5)i 

88 

3.64e-l 

-1.27e-6-(1.22e-5)i 

73 

10 7 

4.16e-2 

-1.32e-7-(1.42e-6)i 

88 

3.94e-l 

-1.34e-7-(1.42e-6)i 

73 


Tables [6791 lists approximation values produced by the CMFE formula Q K ;t[f \, g\ and the FCC formula 
xj. 0 ^ ((/i/</) o g U 1 )) an d computing time they consume. While Table IfTTOl lists approximation values pro¬ 
duced by the CMFE formula Q^° 6 ° 2 [f 2 , g\ and the CFCC formula X^g 1 ^ 12 ((/ 2 / 5 O o gU 1 )) and computing 
time they consume. Both tables show that the CMFE formula consumes significantly less CPU time than 
the FCC, CFCC formulas when they produce comparable approximation results even when the CMFE 
formula uses more functional evaluations. 


7 Concluding remarks 

We develop in this paper composite quadrature formulas for computing highly oscillatory integrals defined 
on a finite interval with both singularities and stationary points. The partitions of the integration interval 
used in the composite quadrature formulas are designed according to the degree of oscillation and the 
singularity. In each of the subintervals, we use piecewise polynomial interpolants to approximate the 
integrand to form two classes of formulas having polynomial (resp. exponential) order of convergence 
by using fixed (resp. variable) number of interpolation nodes. Numerical experiments are carried out to 
confirm the theoretical results on the accuracy of the proposed formulas and to compare them with existing 
methods. Numerical results show that the proposed formulas outperform the existing methods in both 
approximation accuracy and computational efficiency. 
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